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We present results for the corrections of order a 2 (Za) 5 to the Lamb shift. We compute all 
the contributing Feynman diagrams in dimensional regularization and a general covariant gauge 
using a mixture of analytical and numerical methods. We confirm results obtained by other groups 
and improve their precision. Values of the 32 "master integrals" for this and similar problems are 
provided. 

PACS numbers: 31.30.jf, 12.20.Ds 

I. INTRODUCTION II. EVALUATION 



Recent developments in spectroscopy have led to very 
precise experimental values for the 15* Lamb shift and the 
Rydberg constant [l|-[5| , so that now the Lamb shift pro- 
vides the best test of Quantum Electrodynamics for an 
atom. These achievements have spurred great theoreti- 
cal efforts aimed at matching the current experimental 
accuracy (for a review of the present status and recent 
developments in the theory of light hydrogenic atoms, 
see (H). 

The theoretical prediction is expressed in terms of 
three small parameters: Za describing effects due to the 
binding of an electron to a nucleus of atomic number Z; 
a (frequently accompanied by l/7r) from electron selfin- 
ter actions; and the ratio of electron to nucleus masses. 
The Lamb shift is of the order a(Za) 4 ; all corrections 
through the second order in the small parameters are 
known, as well as some of the third order Q. 

Another source of corrections is the spatial distribu- 
tion of the nuclear charge. Even for hydrogen, the ex- 
perimental uncertainty in the measurement of the pro- 
ton root mean square charge radius poses an obstacle for 
further theoretical progress. Fortunately, measurements 
can be performed also with the muonic hydrogen whose 
spectrum is much more sensitive to the proton radius. A 
comparison of the theoretical prediction Q and antici- 
pated new measurements [9( is expected to soon improve 
the knowledge of this crucial parameter. 

In this paper we focus on the second-order radiative- 
nonrecoil contributions to the Lamb shift of order 
a 2 (Za) 5 . The total result for the corrections of this order 
was presented first in [l(| and improved in [lj, [l2| ■ Our 
full result is compatible with the previous ones and has 
better precision. When comparing contributions from in- 
dividual diagrams with [12[, however, we find small dis- 
crepancies in some cases. 

In Section [TT| we present the details of our approach, 
and in Section Hill we present our results. In Appendix lAl 
we show the results for the master integrals, and in Ap- 
pendix [B] new analytic results for two diagrams. 



We consider an electron of mass m orbiting a nucleus 
of mass M and atomic number Z, where Z is assumed 
to be of such a size that Za is a reasonable expansion 
parameter. We are interested in corrections to the Lamb 
shift of order a 2 (Za) 5 and leading order in m/M, given 
by 



5E = -\i) n (0)\ 2 M^' 2 ' 0) (eN -> eN) 



(1) 



where |V>n(0)| 2 = (Za^) 3 / (im 3 ) is the squared modulus 
of the wave function of an S bound state with principal 
quantum number n (fj, is the reduced mass of the system) , 
and M.( 2,2 '°\eN — > eN) is the momentum space repre- 
sentation of the amplitude of the interaction between the 
electron and the nucleus at orders a 2 (Za) 2 and (m/M)°. 
Both particles are considered to be at rest and on their 
mass shell [l3j . 

The correction SE is given by the sum of all the three- 
loop diagrams presented in Figs. [T] and f2J In these fig- 
ures, the continuous line represents the electron, and the 
dashed line represents the interaction with the nucleus. 
The reason for this is that for our purposes this interac- 
tion can be replaced by an effective propagator. In all 
diagrams, the leading order in m/M comes from the re- 
gion where all the loop momenta scale like m. The part 
of the diagrams representing the interaction between the 
electron and the nucleus at order (Za) 2 is given by the 
sum of the direct and crossed two-photon exchange shown 
in Fig. [3j If k and N are the loop and nucleus momenta, 
respectively, and k 2 ~ m 2 <C N 2 — M 2 , the sum of 
the nucleus propagators can be approximated at leading 
order by 
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- M 2 + ie 
iirS(N-k). (2) 



Since the nucleus is considered to be at rest, this gives 
us a S(Mk ). Together with the propagators of the two 
photons, this constitutes the effective propagator. 

We used dimensional regularization, and renormalized 
our results using the on-shell renormalization scheme. 
For all the photon propagators in Figs. [1] and [5] we used 
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Figure 1: The different sets of vacuum polarization diagrams. 
Each set represents the drawn diagram plus all the possible 
permutations of its pieces. 




Figure 3: The sum of the direct and crossed diagrams is ap- 
proximated by an effective propagator (the double line repre- 
sents the propagator of the nucleus) . 
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Figure 2: The diagrams involving a two-loop electron self- 
interaction and vertex corrections. 



a general covariant gauge. The overall cancellation of 
the dependence on the gauge parameter in the final re- 
sult provided us with a good check for our calculations. 
Since the gauge invariance of the sum of the subdiagrams 
in Fig. [3] is trivially and independently fulfilled, for the 
photonic part of the effective propagator we used the 
Feynman gauge. 

Since we are considering free asymptotic states with 
independent spins, the Dirac structures of the electron 
and the nucleus factorize, and we simplify them by in- 
serting each of the structures between the spinors of the 
initial and final states and averaging over the spins of the 
initial states. 

We use the program qgraf [l4| to generate all of the 
diagrams, and the packages q2e and exp [TH, [l6| to ex- 
press them as a series of vertices and propagators that 
can be read by the FORM [ijj package MATAD 3 [l8j|. Fi- 
nally, MATAD 3 is used to represent the diagrams in terms 
of a set of scalar integrals using custom-made routines. In 
this way, we come to represent the amplitude M. in terms 
of about 18000 different scalar integrals. These integrals 
can be expressed in terms of a few master integrals by 



means of integration- by-parts (IBP) identities [19( . Using 
the so-called Laporta algorithm |2| 0| . |2ll| as implemented 
in the Mathematica package FIRE [22], we find 32 master 
integrals. 1 

Since the program FIRE deals only with standard prop- 
agators, when using it we worked only with one of the nu- 
cleus propagators, instead of the Dirac delta of the effec- 
tive propagator. That is, instead of working with §(k ), 
we worked with l/(2k-N+ie), for example. Working with 
just one of the propagators is enough for this purpose, as 
the IBP method is insensitive to the ie prescription (re- 
member from Eq. @ that in our approximation this is 
the only difference between the two propagators). Since 
each diagram in Figs. [T] and [2] represents the subtraction 
of two integrals that only differ in the nucleon propaga- 
tor, when applying the IBP method we can set to zero any 
resulting integral in which the propagator l/(2fc • N + ie) 
disappears. This can be done because the same integral 
with the propagator l/(2k- N — ie) instead would give the 
exact same contribution and thus the difference between 
the two is zero. Once the reduction to master integrals is 
complete, we can simply substitute back the delta func- 
tion in place of the nucleon propagator. 

In order to calculate the master integrals, we turned 
the expressions in Appendix [5] into a representation in 
terms of Feynman parameters. The procedure we then 
followed in most cases was to use a Mellin-Barnes repre- 
sentation [111, H3] to break up sums of Feynman parame- 
ters raised to non-integer powers and transform the inte- 
grals into integrals of Gamma functions over the imagi- 
nary axis. In some cases, we were able to obtain analyt- 
ical results. Otherwise, we used the Mathematica pack- 
ages MB [25[ and MBresolve (2(| to perform a numerical 
calculation. 

For integrals Ig, I w , Iu, hs, hi, hs, hi, and I 32 (cf. 
Appendix [5]) . the Mellin-Barnes representation was too 
cumbersome for a numerical evaluation. In these cases we 
used the Mathematica package FIESTA 1.2.1 [27[ with 
integrators from the CUBA library 28] to perform numer- 
ical computations using sector decomposition [2{| [30| . 
Like FIRE, FIESTA can only process standard propaga- 
tors as input. This means that we had to use the mo- 
mentum representation of the integrals with the nucleon 



1 Actually, it is possible to reduce the number of integrals to at 
least 31, and possibly 30. We give more details about this in 
Appendix lAl 
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propagators instead of the delta function. We separately 
calculated the integrals containing 1/(2N ■ k + ie) and the 
ones containing l/(—2N - k + ie) instead, and added the 
two results. We checked the method by computing with 
FIESTA some integrals we had already found with MB and 
MBresolve. The results always agreed. 

There was one case, integral I32, where the FIESTA 
result for the integral with l/(— 2N - k + ie) was numer- 
ically unstable. Fortunately, in this case we could find 
a representation in terms of Feynman parameters that 
we could compute directly using CUBA, without further 
treatment. This was possible because the integral is fi- 
nite, and the representation was free of spurious diver- 
gences. We cross-checked this result using a beta version 
of FIESTA 2 [3l|, which did not produce the instabilities 
we encountered in the former version. 

We performed an additional cross-check of our results 
by changing the basis of integrals. To do this, we took 
one of the integrals we computed with FIESTA and used 
the IBP method to express it in terms of a similar in- 
tegral of our choice (same as the original one, but with 
some propagator(s) raised to different powers) plus other 
master integrals we already knew. We then computed the 
new integral with FIESTA and checked if the final result 
for the Lamb shift (or for individual diagrams) agreed 
with the calculation in the old basis. Since changing the 
basis modifies the coefficients of all the integrals involved 
in the change, the agreement of the results obtained with 
different bases is a very good cross-check of our calcula- 
tions. 

This cross-check was performed for several integrals. 
In particular, we changed integrals I\g and I27, which 
are the ones limiting our precision, and integrals /15 and 
J32 . Since the last two integrals contain most of the prop- 
agators for integral types F and G, the corresponding 
changes of basis affect the coefficients of most of the other 
integrals of the respective type. 

III. RESULTS 

Our final results for the separate contributions from 
the vacuum-polarization diagrams of Fig. [T]and the dia- 
grams a-s of Fig. [2] are 

a 2 (Za) 5 / u \ 3 
SE vac . = { ™ (£■) m [0.86281422(3)], (3) 

6E a -„ = ?l¥*L (JL) 3 m [-7.72381(4)] . (4) 

The best results so far for the vacuum polarization dia- 
grams and for diagrams a-s have been published in [32| 
(cf. [1] for references of partial results) and [l2[ , respec- 
tively. Our results are compatible with them and improve 
the precision by two orders of magnitude in the case of 
$E vac . and a little over one order of magnitude for 5E a - s . 
The total result reads 

SE = a ^^ Z( f (V_) 3 m [-6.86100(4)] , (5) 



Table I: Comparison between our results for the different 
vacuum-polarization sets (in Fried- Yennie gauge) and those 
of [13, • Numbers ending in an ellipsis indicate an analytic 
result, which we show in Appendix iBl 

Set This paper Refs. [32, 35] 

I -0.07290996446926(4) -0.0729098(3) 

II 0.61133839226... 0.61133839226... 

III 0.50814858506. . . 0.50814858506 . . . 

IV -0.12291623(3) -0.122915(3) 

V -23/378 -23/378 



Table II: Comparison between our results for diagrams a-s 
(in Fried- Yennie gauge) and those of [l2T ]. Numbers ending 
in an ellipsis indicate an analytic result, which we show in 
Appendix [B] 



Diagram This paper Ref. [12] 

a 

b 2.955090809... 2.9551(1) 

c -2.22312657... -2.2231(1) 

d -5.2381153272259(2) -5.238023(56) 

e 5.0561650638185(4) 5.056278(81) 

/ 61n2 - 207/40 -1.016145(21) 

g 6 In 2 - 147/80 - tt 2 /4 -0.1460233(52) 

h 153/80 153/80 

i -5.51731(2) -5.51683(34) 

j -7.76838(1) -7.76815(17) 

k 1.9597582447795(2) 1.959589(33) 

I 1.74834(4) 1.74815(38) 

m 1.87510512(6) 1.87540(17) 

n -1.30570289(7) -1.30584(18) 

o -12.06904(9) -12.06751(47) 

p 6.13815(1) 6.13776(25) 

q -7.52425(2) -7.52453(34) 

r 14.36962(7) 14.36733(44) 

s -0.9304766935602(5) -0.930268(72) 



and the corresponding energy shifts for the IS and the 
2S states in hydrogen are 

SE 1S = -296.866(2) kHz, (6) 
5E 2S = -37.1082(3) kHz. (7) 

Choosing the Fried- Yennie gauge [33|, HH , we also com- 
pared the results from the different sets of vacuum po- 
larization diagrams with those of [32|, [3|| , and the results 
from the individual diagrams a-s with those of [l2} . Our 
results for the vacuum polarization graphs and diagrams 
a-s are presented in Tables U and HH respectively. All 
numbers in the tables are to be multiplied by the prefac- 
tor a 2 (Za) 5 / (7ro 3 )(/x/m) 3 m (note the difference in nor- 
malization in [32j). 

We found new analytic results for four diagrams. The 
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results for diagrams / and g are given in Table HH while 
the results for diagrams b and c, being too lengthy for 
the table, are presented in Appendix |E1 For complete- 
ness, the known analytic results for sets II and III of the 
vacuum polarization diagrams are given in Appendix [B] 
as well. 

It should be mentioned that the errors of the results 
in Eqs. ©, (jU), and (JSJ are not obtained from the sum 
of the errors of the diagrams in Tables Q] and [TTJ Once 
we decompose the problem into the calculation of mas- 
ter integrals, the diagrams are no longer independent, as 
the same master integral contributes to several different 
diagrams. Thus, to find the error of our total result, we 
first sum all diagrams and then sum all the errors of the 
integrals in quadrature. 

We found discrepancies between our results for dia- 
grams as and those of [l2j. Most of the central values 
in the second and third column of table|TT]lie between la 
and 2(7 away from each other, but in the case of diagrams 
o and s the difference is around 3c, and for diagrams k 
and r, it reaches 5a (we take as a the errors of individual 
diagrams in the third column). We should stress again 
that our calculation is done using dimensional regular- 
ization while the study of Ref. [l2| was performed in four 
dimensions. Even though all the individual diagrams are 
finite, one can imagine situations where the two regular- 
ization methods give different partial results. However, 
we do not observe significant cancellations in the sum of 
the differences. Thus, it seems the differences are real 
although practically negligible; their sum is very small 
and amounts to 10 -3 , which is the error estimate in fl2l . 



Thus our results agree within that error. 



IV. SUMMARY 

We have applied particle theory methods to compute, 
in dimensional regularization and a general covariant 
gauge, the corrections of order a 2 (Za) 5 to the Lamb 
shift. We have made use of IBP techniques to reduce 
the problem of computing all the necessary Feynman di- 
agrams to the simpler problem of computing 32 scalar 
integrals. Mellin-Barnes integral representations and sec- 
tor decomposition have then allowed us to obtain analytic 
results for some of these integrals, and good numerical re- 
sults for the rest. With this, we have been able to repro- 
duce and improve the results from previous calculations. 
The techniques used here are quite general and can be 
applied to other multi-loop problems in atomic physics. 
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Appendix A: Results for the master integrals 



In section |TT] we presented our method of calculating the corrections to the Lamb shift, which differs significantly 
from the methods used in previous calculations. One important difference is the reduction of diagrams to master 
integrals. Here we present our results for all master integrals. 

The set of master integrals is represented in Fig. |H We have two different types of integrals. In Euclidean space, 
they are defined as: 



e 37E« r 

F(l>i,l>2,l> 3 , 2/ 4 , 2/5, U 6 , V 7 , V%) = 3 / 

(tt d / 2 ) J 



d D h d D k 2 d D k 3 27T(5(fc§) 

(fc2)"i (fc2)"2 (jfc2)„ 3 +p y + ip; 

1 



[(fci + k 2 + p) 2 + I]"* [(fci + k 2 + k 3 + p) 2 + 1]* [(k 2 + k 3 + p) 2 + [(k 3 + pY 

e 3 ^ e f d D k 1 d D k 2 d D k 3 2ir5(k<l) 



G{vi,v 2 ,v 3 ,v i ,u b ,v 6 ,v 7 ) 



J (fcf)"i (k\y-> (k$) u » [(fci + k 2 +pY + l}" 4 

1 



[(fci + k 2 + k 3 + p) 2 + l]" 5 [(k 2 + k 3 + pY + 1]^ [(fc 3 + PY + I]" 7 ' 



(Al) 



(A2) 



where D — 4 — 2e, and p — (i, 0) is the momentum of the electron. The mass of the electron has been set equal to 
one for convenience and can be easily restored from the dimension of the integral. The factor e 37E£ , where je is the 
Euler-Mascheroni constant, has been introduced to suppress the dependence of the results on this constant. 
With these definitions, our results for the master integrals are: 

h = F(l, 0, 0, 0, 0, 1, 0, 1) = 2e 3 ^ -i V 2 r( _ 3 ^ 4e) " > (A3) 
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Iqs he hi hs ho ho hi 



h2~ 

Figure 4: A graphic representation of the 32 master integrals. Solid and dashed lines represent massive and massless scalar 
propagators, respectively. The dotted double line denotes the delta function. A dot on a line signifies that the propagator is 
raised to a higher power. The external lines indicate the momentum p of the electron flowing in and out of the diagram. The 
first two diagrams represent different integrals that differ only by a term in the numerator. 
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F(l,0,-l,0,0,l,0,l) = -2Ji, 

F(0, 0, 0, 0, 1, 1, 1, 0) = -96.174642407742494299(1) - 3003.97283051743374945(l)e 

-16370.644886761701890(l)e 2 - 204040.09217878970569(l)e 3 + C(e 4 ) , 

F(-l, 0, 0, 1, 0, 1, 0, 1) = 128.23285654365665907(1) + 4005.2971073565783326(l)e 

+21827.526515682269186(l)e 2 + 272053.45623838627426(l)e 3 + 0(e 4 ) , 

F(0, -1, 0, 0, 1, 1, 1, 0) = 213.37528929773859515(1) - 1789.0076495990746772(l)e + 0(e 2 ) , 



F(l, 0,1, 0,0, 1,0,0) = 2^Fe 
F(l, 0,0, 0,1, 0,0,1) = 2^e 



: „,„ r 2 (i-e)r(-f + 3e)r(-f + 2e)r(f-5e) 

r(2-2e)r(3-3e) 



•2e)r(|-3e)r(-i + e ) 



F(0, 0,0, 0,1, 0,1,1) = 2V?re 

8^ 2 



3 ^ 6 r(-i + £ )ri 

r 

r(-i + e)r(-f + 2e)r 2 (-| + e) 



T(2 - 2e) 



F(0, 1,0, 0,1, 1,1,1) 
F(0, 0,0, 1,1, 1,1,1) 

F(l, 0, 0, 0, 1, 1, 0, 1) = 2 3 - 4e 7re 37Ee 
_r(l - e)T (3e 



r(-l + 2e) 

— - 257.35053226188(1) - 2952.9668342496406(4)e + 0(e 2 
420.49901(1) + 1860.837(4)e + 0(e 2 ) , 



cos (27re) 



r(|-3e)r(e) / 5 3 

2 _V2^ >_±1 p 2 1, - -3e,e;-,4-4e;l 

r (4 - 4e) V 2 2 ' 



3 5 



— m s — T-r~ 3F2 1, 1 - e,3e - -; - - 2e, 2e; 1 

r(|-2e)r(2e) 3 V 2' 2 



4tt 



An 
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F(l, 1,0, 0,1, 1,0,1) = — - 24tt 2 - — +7T 2 ^-116- yTr 2 + 321n2 + 100C(3)J e + C(e 2 ) 

F(l, 1, 0, 1, 0, 1, 0, 1) = -263.74028719521945979(1) + 1741.1125810306205720(l)e + 0(e 2 ) , 
F(l, 0, 0, 0, 1, 1, 1, 1) = -362.8560(1) + 0(e) , 
F(l, 1, 0, 0, 1, 1, 1, 1) = 36.969282(2) + 0(e) , 



(A4) 

(A5) 

(A6) 
(A7) 

(A8) 

(A9) 

(A10) 

(All) 
(A12) 

(A13) 

(A14) 

(A15) 
(A16) 
(A17) 
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he = F(l,4,0,0,l,0,l,l) = 7r J - 



'2.-, 



1 27 

h% 
h<) 



hi 

hi 



2I 513 H077 „ 571 ^^^/iWS 



hi = 
hs = 

h9 = 

ho = 
hi = 



128e 768 16 V 2 

1889.3810189605726842(l)e - 2199.2559561980712031(l)e 2 + 0(e 3 ) , 

683.43054120051764110(l)e 



F(l, 0,0, 0,2, 0,2,1)=^ In' ' X " 



+5647.2496334930969112(l)e 2 + 0(e 3 ) , 
F(1,4,1,0,1,0,1,1) = tt 2 
F(l, 0, 1, 0, 1, 1, 1, 0) = -293.4480(2) + 0(e) , 
F(l,0,0,l,l,0,l,l) = 7r 2 



343 125257 2 9 169 , „ „, 2 „ , / l + \/5 

1 7r H n 2 + 16 n 2 2 - 48 n 2 — 

512e 46080 3 192 \ 2 



-2-2 
e 



32 In 



1 + V5 N 



+1394.0754186124348755(l)e+ C(e 2 ) , 

47T 4 

F(l, 1, 1, 0, 0, 1, 0, 1) = 2tt 2 - — + ir 2 (44 - 4tt 2 + 80C(3)) e + 0(e 2 ) , 



hi = F(0, 0,0, 1,1, 0,1,1) = - 
/ 23 



1287T' 



F(l, 0,1, 0,1, 0,1,0) 

27r 3 e 3 T Ee 
T 2 (2-2e)r(| -e) 



V 3 

2e _ 2 r(i-e)r( e -i 



1792 256 1024, , . , 
■ — n-— In 2 e + 0(e 2 ), 



(A18) 

(A19) 
0(e) , (A20) 
(A21) 

(A22) 
(A23) 

(A24) 



r(|-5e) sin (ire) 



T (3e- 



sin (ire) cos (47re) \T (§ - 3e) sin (2ire) cos (37re) r (5e - |) cos (2ire 



V^r(2-2e)r(2e- §) 



r (2 - e) T (3e - |) sin (ttc) cos (?re) cos (3?re) 



F 2 I l,2-2e,2e- i;2-e,3e- 



2 1 - 2 ^r(|-3e) 



2e) T (i + e) sin (2ire) cos (ire) cos (27re) 



iF, ( l,^-3e,e;^-2e,2e;l 



hi = G(0, 1,2, 1,0, 1,0) 



2tt 2 



162.745878930257(1) + 640.681562239(2)e 



-9490.745115169417(3)e 2 + 0(e 3 ) 
Air 2 



G(2, 1,1, 1,0, 1,0) 



192.3546921335253(1) - 2297.18352848038(l)e 



-10356. 58582995624(l)e 2 + 0(e 3 ) , 
4^ 2 



G(l, 1,1, 1,0, 1,0) 



- 244.4995291143211(3)- 2339.54007847666(2)e +Q(e A ) , 



G(0, 1, 1, 2, 0, 1, 1) = 136.8086023(2) - 907.048(2)e + 0(e 2 ) , 
G(0, 1, 1, 1, 1, 1, 0) = -280.62418(1) + 734.494(l)e + 0(e 2 ) , 
G(l, 1, 1, 0, 1, 1, 1) = 118.63826101784(1) + 0(e) , 



ho = G(0,1,1,0,1,1,0) = - 



V r(2-2e)r(|-4e) 



3e 



X3 J F 2 Q-5e,^-e,-i + e;^-4e,i + e;l 

G(l, 1, 1, 1, 1, 1, 0) = 49.3616(1) + 0(e) , 

G(l, 1, 1, 1, 1, 1, 1) = 26.272804(6) + 291.1097(l)e + 0(e 2 ) , 

I 



(A25) 

(A26) 

(A27) 
(A28) 

(A29) 
(A30) 
(A31) 

(A32) 

(A33) 
(A34) 



where £ denotes Riemann's zeta function, and 3F2 is a 
generalized hypergeometric function. The latter can be 
expanded in e with the help of the Mathematica package 



HypExp 2 [38j. 

The relation between integrals h and h expressed in 
Eq. (IA4[) is not evident when looking at their respec- 
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tive representations. This relation becomes clear when 
checking the cancellation of the gauge-parameter depen- 
dence in the sum of all diagrams. If the 32 integrals 
presented here were an irreducible basis, the gauge de- 
pendence of the coefficient of each integral should vanish 
independently. However, this does not happen with the 
coefficients of integrals I\ and h, which means that the 
integrals are connected. Demanding the cancellation of 
the gauge dependence yields Eq. (IA4|) . We checked this 
relation by computing explicitly the analytic solution for 

There appears to be also a relation between integrals 
I3 and I4, although the gauge dependence does not give 
us any hint in this case. By demanding the cancellation 
of poles in several diagrams, one can find the following 
relation between the first three terms of J3 and I4, 



tioned above, we were able to numerically calculate these 
integrals to 100-digit precision and confirm the analytic 
expressions with PSLQ. 



Appendix B: Analytic results 

Here we show the analytic results for diagrams b and 
c from Fig. [2j 



Diagram b = 
111 



tx a -91n2 + 241n 2 



5- fa*. 



-72 In 



1 + V5 



(Bl) 



3 



h + 0(e 3 ) . 



(A35) 



The relation, however, seems to be valid to all orders in 
the e expansion. We checked it numerically up to order 
e 4 , but we could not find an analytic proof for it. 

Integrals I3-I5, I12, I13, he, hi, and ho~hi can be 
represented as a one-fold Mellin-Barnes integral. We only 
show numerical results with 20-digit precision, which is 
more than enough for our purposes. However, these in- 
tegrals can be easily evaluated with a precision of 100 
digits or more. With this kind of precision it is possible 
to find analytical results, using the PSLQ algorithm [39| . 
In this way, we determined the 0(e) term of integrals I12 
and hi- 

The analytic expression for the O(e ) term in ho was 
obtained using the analytic result for set II of vacuum- 
polarization diagrams presented in [35[ . Likewise, the an- 
alytic expression for the 0(e) term in I22 was extracted 
from the analytic result for set III found in [32| . As men- 



Diagram c = 

352897 31 , 643 , n 248 , 2 „ 

1 7T In 2 hr 2 

27000 45 225 15 



104 , (l + y/5 
: In 



248 



In 2 



1 + V5 



(B2) 



9V5 "\ 2 J 5 "' \ 2 / 
For completeness, we also give here the analytic results 
for sets II and III of the vacuum polarization diagrams, 
found in [35j and [32j . respectively: 



Set II = 



67282 2 
6615 ~ 9' 



628 
~63~ 



In 2 



Set III 



-^V51n(i±^ 
63 \ 2 



15647 25 52, „ 

ttH ln2. 

13230 63 63 



8 ln2 /lW5 
3 I 2 
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